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' We consider the statics and dynamics of F = 1 spinor Bose-Einstein condensates (BECs) 

confined in double well potentials. We use a two-mode Galerkin-type quasi-analytical approx- 
imation to describe the stationary states of the system. This way, we are able to obtain not 
only earlier results based on the single mode approximation (SMA) frequently used in studies of 
spinor BECs, but also additional modes that involve either two or all three spinor components 
of the F — 1 spinor BEC. The results based on this Galerkin-type decomposition are in good 
■ agreement with the analysis of the full system. We subsequently analyze the stability of these 

^ . multi-component states, as well as their dynamics when we find them to be unstable. The 

instabilities of the symmetric or anti-symmetric states exhibit symmetry-breaking and recur- 
rent asymmetric patterns. Our results yield qualitatively similar bifurcation diagrams both for 
polar (such as ^"^Na) and ferromagnetic (such as ^^Rb) spinor BECs. 
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Abstract 



1 Introduction 

In the past few years, there has been a remarkable amount of progress in the study of Bose-Einstein 
condensates (BECs) [U [2] and an intense investigation of the localized nonlinear states that can be 
formed therein [3j . Prom the nonlinear dynamics point of view, one of the features that makes this 
system a particularly relevant and interesting one is associated with the presence of a diverse range 
of external potentials which are used to trap the atoms magnetically, optically or electrically (or 
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through combinations thereof) [HEllS]. Hence, monitoring the existence, stabihty and dynamical 
properties of nonhnear locahzed modes within these diverse potentials has become a principal 
theme of research effort, especially within the realm of the prototypical mean- field model, namely 
the Gross-Pitaevskii (GP) equation (which is a variant of the nonlinear Schrodinger (NLS) equation 
extensively used in nonlinear optics [4]). 

While most of the work in BECs has centered around one-component systems, more recently far- 
off-resonant optical techniques for trapping ultracold atomic gases [5], have produced an intense 
focus on the study of spinor BECs |6l [7] , in which the spin degree of freedom (frozen in magnetic 
traps) emerges. In this context, various phenomena absent in single-component BECs may arise 
in the dynamical evolution of multi-component spinor condensates [8], including the formation of 
spin domains [9] and spin textures [10], spin-mixing dynamics [TT], dynamic fragmentation [12j, or 
dynamics of quantum phases [13]. Moreover, macroscopic nonlinear states in the form of multi- 
component vector solitons of the bright [HI [151 [16] , dark [17] and gap [18] type, along with more 
complex structures, such as bright-dark soliton complexes [19] and domain walls [20], have also 
been studied. 

On the other hand, as indicated above, one of the most attractive traits of such ultracold systems 
is the possibility of experimental realization of different potentials. Among them, one that has 
drawn considerable attention due to its fundamental nature is the double-well potential. One of 
the prototypical realizations thereof comes from combining a strong harmonic trap with a periodic 
lattice |2j. In this context, the study of BECs loaded in double well potentials allows for the 
investigation of a variety of fundamental phenomena, including Josephson oscillations and tunneling 
for a small number of atoms, or macroscopic quantum self- trapping and an asymmetric partition 
of the atoms between the wells for sufficiently large numbers of atoms [21j. Double well potentials 
have also spurred numerous theoretical insights including, among others, finite-mode reductions, 
analytical results for specially designed shapes of the potential and quantum depletion effects [22l 
[23[Mll25l|26l[27l[28l[29l|30^ It should be noted that such potentials have also been studied 
in the context of nonlinear optics, with relevant experimental results appearing, e.g., in twin-core 
self-guided laser beams in Kerr media [31] and optically induced dual-core waveguiding structures 
in photorefractive crystals 

The principal scope of the present work is to combine these two interesting developments, namely to 
examine spinor BECs in the presence of double well potentials. In particular, our aim is to present 
a systematic classification of the states that are possible for F = 1 spinor condensates confined 
in a double well potential. This classification is formulated on the basis of a two-mode Galerkin- 
type approximation of the stationary states of the system, which involves the decomposition of 
the spatial part of the solutions into a linear combination of the eigenfunctions of the underlying 
linear operator; the relevant two- mode reduction used here is in the spirit of Galerkin truncations 
in finite element and similar methods (see, e.g., Ref. [33j). Such an approach can provide a detailed 
analytical handle on the dynamics of this multi-component BEG system (see, e.g., the earlier work 
on one-component [29j and two-component BECs [M] in double well potentials). It should be 
noted here that such attempts have been made before, most notably in Refs. [35] and [36]: the 
first of these works examined magnetization oscillations and beats among other dynamical states, 
while the second one (which provided a description beyond mean-field theory) examined quantum 
entanglement and pseudo-spin-squeezing properties. However, these works were constrained within 
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the commonly used single mode approximation (SMA) for spinor condensates (see, e.g., p!T| [37| [38] ) 
in which the wavefunction of each of the three components is taken to be a time-dependent multiple 
of a single stationary state. One of the particularly interesting features in the present setting is 
that the two- mode Galerkin-type approximation reveals a considerable wealth of states, many of 
which can not be described in the framework of SMA. In particular, as we will see below, there is a 
multitude of states which are symmetric in some of the components, while they are anti-symmetric 
in others, and out of these states also bifurcate further states which are asymmetric (but often in a 
non-single- mode way) in the different components. Many of these novel states have a good chance to 
be observed in experiments with quasi one-dimensional (ID) spinor BECs since, by performing their 
linear stability analysis, we find them to be stable. Finally, we examine the dynamical evolution 
of both two-component and three-component states that we find to be unstable; we observe that 
the manifestation of the associated instabilities is exhibited typically through symmetry-breaking 
between the components and is also associated with recurring asymmetric patterns. 

The paper is structured as follows. In section 2, we present the model, and provide the analytical 
approach. In section 3, we present our numerical results. In particular, we obtain the complete 
bifurcation diagram of the possible stationary states, both for the full three-component system, 
as well as for the two-mode Galerkin-type approximation (in very good agreement between the 
two). In that section, we also illustrate the spatial profiles of the nonlinear modes arising in the 
spinor condensate and quantify their stability. Lastly, we corroborate these stability predictions by 
performing dynamical simulations of the unstable modes. Finally, in section 4, we summarize our 
findings and present our conclusions, as well as some topics for future study. 



2 The model and the analytical approach 

In the framework of mean-field theory, the wavefunctions ip±i,o{x^t) of the three hyperfine compo- 
nents {rriF = il,0) of a quasi-lD spinor F = 1 condensate is governed by the following system of 
coupled normalized GP equations [18]: 

idt^±i = C^±i ^ iysS^±i ^ MS - 2\^^i\^)^^i ^ ua^lr^, (2.1) 

In these expressions S = IV^-iP + IV^oP + IV^iP represents the total normalized density, while the 
coupling coefficients and Ua represent, respectively, the symmetric spin-independent and the 
antisymmetric spin-dependent interaction strengths (see, e.g., [18j for details). Additionally, 

C^j = -\d^^^j+V{x)4>j, (2.3) 

{j = —1, 0, 1) represents the single-particle operator with a confining potential V{x) assumed to be 
of the form 

V{x) = ^n^x^ + Vosech^{x/w). (2.4) 

This is a double well potential consisting of a parabolic potential of strength Q emulating, e.g., a 
usual harmonic trap, and a localized barrier potential of strength Vb and width representing, e.g., 
a blue-detuned laser beam that repels atoms from the harmonic trap center. Notice that the main 
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qualitative features of our results (i.e., the nature of the bifurcation diagram near the linear limit, 
the emergence of the symmetry-breaking asymmetric states, and their dynamical manifestations) 
do not rely on the specific form of the double well potential and should be expected to arise more 
broadly within this class of models (spinor condensates in double well potentials). 

We now seek stationary solutions of Eqs. (|2.ip - (|2.2p in the form = uj ex.p{—ifijt) exp(i6>j), where 
jiij and Oj correspond, respectively, to the chemical potentials and phases of the wavefunctions 
obeying the usual constraints (see, e.g., [37l[20]), namely, 2/io = Mi+M-i and AO = 20o — {0i-\-0-i) = 
or TT. These are the so-called phase matching conditions relevant for stationary states typically 
in systems with parametric interactions; see, e.g., Ref. [39] for the case of four- wave mixing in 
nonlinear optics. This ansatz results in the following equations for the stationary states Uj: 



(2.5) 
(2.6) 



where the ± sign in the last terms of the above equations corresponds, respectively, to the cases of 
A0 = and AO = tt. 

The basic idea of the two mode Galerkin-type approximation [23l [26l [28l |29] that we will employ 
to approximate the solutions of Eqs. (|2.5p - (|2.6p is the following: we assume that in the vicinity of 
the linear limit of the system (and for appropriate selection of the chemical potential) each of the 
states u±i^o can be decomposed into a two- mode expansion, involving the symmetric ground state 
and the antisymmetric first excited state of the underlying linear potential. We will denote those 
states by (j)a and 05, respectively, and accordingly decompose Uj (where j G { — 1,0, 1}) as follows. 



(2.7) 



where Ca^ and c^-^^ are unknown time-dependent complex prefactors. By substituting Eq. (|2.7p into 
Eqs. (|2.5p - (|2.6p . and upon projecting each of the three equations in (|2.5p - (|2.6p to (j)a and 05, we 
obtain through tedious but straightforward algebra, the following six nonlinear, algebraic equations 
describing the stationary states of the original system 



(±1) 

/^±1C5 



± 



± 



(±1) 



(±1) 



.(±1)' 



(±1) 



(2.8) 



C (±1) 

^bSbcl 



(2.9) 



(0) 



TaiSa 



-"^Note that since we are interested in stationary states herein, we will consider and cj^^ to be independent of 
time, although the expansion can also be used to probe the dynamics, by substituting the ansatz of Eq. (j2.7|) with 
time-dependent coefficients into Eqs. (|2.1|) - (|2.2|) and projecting into (f)a and (f)i,. 
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± 



(2.10) 



,(0) 



± 



(2.11) 



In these equations, uja^h are the eigenvalues that correspond to the hnear eigenfunctions ^^,6, re- 



spectively, Sa = E,(ci'0', ^6 = E,(4 = 2ci^'^4^\ ^hilc the coefficients = / <i>\dx, 



= / ^^(ix and r^b = / (\)\(\)\dx are constants that depend on the potential. 

By solving this algebraic set of equations, for appropriate parameters of the potential [e.g., in the 
normalization detailed below (see Section 3) using the mean- field interaction energy, for ^ = 0.2, 
Vo = l Sindw = 0.5, Ua 0.249, 0.315, and Ta 0.127, 0.134 and Tab 0.120], we can 
extract all the potential stationary states for a given combination of chemical potentials in the two- 
parameter space (/io^Mi) [recall that is fully determined by the above two parameters]. Then, 
these can be used as initial guesses for identifying these solutions in the full nonlinear eigenvalue 
problem of Eqs. (1231) - (12:61) . 

Some special cases of solutions can be found in the above algebraic equations. Since these have been 
discussed in earlier works in one-component systems (see, e.g., [29]), two-component systems (see, 
e.g., [S]), and even in the full three-component system but in the context of the SMA (see e.g., 
[m [m [33 [38]), we relegate the relevant detailed analysis of these to Appendix A. Importantly, 
we note that the principal building blocks of the single component setting are a symmetric state, 
where only Ca (of a particular component) is non-zero, an anti-symmetric state with only C5 7^ 0, 
as well as asymmetric states with both Ca and C5 7^ 0. The latter bifurcate beyond a critical point 
from the symmetric state if < 0, or from the antisymmetric state if > [29] , 

In the next section, we will offer the full set of solutions that can be constructed within our two-mode 
Galerkin-type approximation and compare them to the detailed numerical results of the original 
problem. 

3 Numerical Results 

Before proceeding further, it is necessary to provide here values of the physical parameters that we 
will use. First we note that for simplicity in our computations we use the rescaling of Vs to 1 and of 
Ua to i^a/^s = Notice that in the relevant cases of the ferromagnetic ^^Rb and polar ^^Na spinor 
{F = 1) condensates, the parameter S takes the values S = —4.66 x 10~^ [40] and S = +3.14 x 10~^ 
[31], respectively. Under this rescaling, the connection between physical values and dimensionless 
parameters is as described in detail in Refs. [191 [SO]- We will also fix the double well potential 
parameters to the values 1] = 0.2, = 1 and w = 0.5 (other values lead to qualitatively similar 
results). Taking into regard that the normalized chemical potentials will take values /ij ~ 0.5, this 
choice may correspond, e.g., to a spinor condensate of sodium atoms confined in a highly asymmetric 
trap with frequencies uj± = 3uJx = 27r x 230 Hz, with peak ID density no — 10^ m~^ and number 
of atoms of order O(IO^). In this case, the time and space units in the numerical results that will 
be presented below are 1.2 ms and 1.8 /im, respectively. 
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Figure 1: (Color online) Top panels: the norm N (normalized number of atoms) of the numerically 
found solutions of Eqs. (|2.ip - (|2.2p (left) and their counterparts predicted by the two- mode Galerkin- 
type approximation (right) for the case of the ^^Na spinor BEC, as a function of /io with fixed 
= 0.38. Middle and bottom panels: blowups of segments in the top left panel where new 
combined solutions emerge at the relevant bifurcation points. All blowups, except the last one, are 
dedicated to the case of = 0.38 as in the top left panel. The last one, showing the branches C7 
and C8 (see text) corresponds to the case of = 0.48 (as in the top left panel of Fig. [71 below). 
The branches Dl and D3 are generated from D2 and D4, respectively. Branches labeled by D (Dl- 
D6 in the different panels of the figure) correspond to waveforms with two nonzero components, ui 
and U-i; see the discussion of section [311 and Figs. [3]-[l]for more details. The branches CI, C3, C5 
and C7 are generated from C2, C4, C6 and C8, respectively. Branches labeled by C (C1-C8 in the 
different panels of the figure) correspond to waveforms with all three components nonzero; see the 
discussion of section 13.21 and Figs. [5]l6] for more details. Solid (blue) lines and dashed (red) lines 
denote stable and unstable solutions, respectively. 
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Perhaps the cornerstone of the present work is encapsulated within Fig. [TJ generated for the case 
of a ^^Na spinor condensate. This, admittedly, rather busy diagram encompasses all possible states 
near the linear limit of a spinor BEC loaded in a double well potential. Among them, one can discern 
one-component states that are well-known from the considerations of single-component double well 
settings [21] [22l |26l [29] for each one of the three components (labeled by S for symmetric, AN for 
antisymmetric and AS for asymmetric realizations of each component); see also Fig. [2l They are 
labeled by 1 when they belong to component iii , by when they belong to i^o and by 2 when they 
belong to U-i. Notice that the AS state bifurcates from the S state in the attractive case and from 
the AN state in the repulsive case, inheriting their stability (and rendering them unstable, through 
a pitchfork bifurcation) [29j. 

Importantly, the graph also contains two-component states involving components with subscripts 1 
and —1 (and with a vanishing component), denoted collectively by D and analyzed in Section \3A\ 
below (cf. also [SJI)- These are branches involving the symmetric and/or antisymmetric components 
of Ui and U-i; interestingly, out of these branches bifurcate new asymmetric solutions involving the 
same two components. 

Finally, branches with contributions from all three components are also identified and are collectively 
labeled by C; these are examined more systematically in Section [3^ below. These predominantly 
consist of either components ui and uq (with a small component inu-i) or U-i and uq (with a small 
component Ui). Again different combinations of symmetric and/or antisymmetric configurations 
are possible among these components, and additional asymmetric branches are observed to bifurcate 
from them. We now turn to a more detailed explanation of our findings. 

To analyze the stability of a particular state (in this case and in the following ones), we perform 
linearization around the unperturbed state Uj^ assuming a perturbed solution of Eqs. (|2.ip - (|2.2p in 
the form, 

jpj = exp{-ijj.jt) [uj + e {pj{x) exp(At) + qj{x) ex.p{X^t))] (3.1) 

where pj and qj represent infinitesimal perturbations with eigenvalues A = Xr -\- iXi (here e is a 
formal small parameter and the asterisk denotes complex conjugate). When its relevant eigenvalues 
are found to be purely imaginary, then the corresponding state is linearly stable, while the presence 
of eigenvalues with 7^ is tantamount to instability. Then, the results of the linear stability 
analysis are typically presented in the spectral plane (A^, A^). A typical example is provided in Fig. 
[21 where we show the profiles and the corresponding spectral planes (Ar,Ai) of the branches SI, 
ANl, and ASl. 

Coming back to the bifurcation diagram of Fig. [H it is important to highlight that it has been 
constructed by fixing and varying only /ii and /io (since one of these three quantities is always 
slaved to the variation of the other two). For this reason, S2, AN2 and AS2 are horizontal lines 
since is fixed in the case studied below (i.e., their properties do not change as (/ii,/io) are 
varied.) We will vary /io while doing the analysis, with /ii being determined by fii = 2/io — M-i- 

3.1 Two- component states 

In addition to "pure" branches populating only one component, there also exist additional branches 
populating two components, similarly to the results of Ref. [34]. In Fig. [U prototypical examples of 
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Figure 2: (Color online) Profiles of wave functions ui (first and third rows) and stability eigenvalues 
(second and fourth rows) corresponding to SI, ANl and ASl (branches of these solutions are shown 
in Fig. [1]): SI for /io = 0.48 (top left), ANl for /io = 0.48 (top right), ASl for fio = 0.48 (third row 
left) and ASl for /io = 0.45 (third row right). 



these branches are given by D1-D6 which involve the two components Ui and U-i. More specifically, 
D2 "connects" SI and AN2. That is to say that D2 has a symmetric wavefunction in the component 
iii, while it has an antisymmetric one in component U-i. Dl bifurcates from D2, below a critical 
number of atoms through a pitchfork bifurcation, is asymmetric in both components ui and U-i 
and eventually merges into AS2. It is relevant to point out here that in the first panel of the second 
row of Fig. [U it appears as if D2 and Dl are essentially overlapping for the small fraction of the D2 
branch for which they coexist (in this segment, D2 is unstable). However D2 terminates into AN2, 
while Dl continues to exist for lower norms, down to AS2. In a situation entirely similar to that 
highlighted above, but now between S2 and ANl, the branch that "connects" them is D4. Out of 
that branch with a symmetric wavefunction in U-i and an anti-symmetric in ui bifurcates (above 
a critical norm A^) the branch D3 which is again asymmetric in both of these components, and will 
eventually merge with ASl. Since these two situations (involving D1,D2 and D3,D4, respectively) 
are essentially similar (both of the bifurcation sub-structures are shown in Fig. [T]), only the latter 
pair of branches is examined in detail in Fig. [3l In particular, it can be seen there that while the 
D4 branch is stable before D3 appears, the supercritical pitchfork leading to the emergence of D3 
destabilizes D4. On the other hand, D3 itself may be unstable but only due to Hamiltonian-Hopf 
bifurcations [32], associated with weak complex quartets of eigenvalues. Finally, in terms of two- 
component branches, in addition to the above pairs connecting the symmetric nonlinear eigenmodes 
of one component with the anti-symmetric ones of the other, there also exist the branches D5 and 
D6 (see again the blowups of Fig. [1]). These connect the symmetric state of ui (branch SI) with 
the symmetric state of U-i (branch S2) and the antisymmetric state of ui (branch ANl) with the 
antisymmetric state of U-i (branch AN2). The branch D5 is linearly stable, while D6 is linearly 
unstable, as indicated in Fig. [H Interestingly, the norm of these branches almost coincides with 
the norm of branches SO (for D5) and ANO (for D6) although they do not involve a nontrivial 
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component of uq. This is apparently associated with the fact that branches such as D5 or D6 can 
be effectively described in the framework of SMA. On the other hand, clearly D1-D4 cannot; thus, 
they provide a distinct, potentially stable (and even symmetry-broken) set of states that could be 
accessible by the F = 1 spinor BEG system. 




Figure 3: (Color online) Left panels: the maximum real eigenvalue as a function of fio for D3 (top 
left) and D4 (bottom left) in Fig. [H The middle and right panels show the profiles of wave functions 
(first and third row) ui and U-i [by solid (blue) and dash-dotted (red) lines, respectively], together 
with the stability eigenvalues corresponding to D3 and D4: D3 for /io = 0-48 (top middle), D3 for 
/io = 0-51 (top right), D4 for /io = 0-41 (bottom middle) and D4 for /io = 0.42 (bottom right) with 
/i_i = 0.38 in the case of ^^Na condensate. 




Figure 4: (Color online) The left panel shows the maximum real eigenvalue as a function of /io for 
D6 in Fig. [TJ The middle and right panels show the profiles of wave functions ui and U-i by solid 
(blue) and dash-dotted (red) lines, respectively, together with stability eigenvalues corresponding 
to D5 and D6 in Fig. [B D5 for /io = 0.382 (middle), D6 for /io = 0.381 (right). 



3.2 Three-component states 

We now turn to states involving all three components, such as the ones analyzed in Figs. [5]and[6l 
These branches resemble the ones that we analyzed above regarding states D1-D4; in particular, 
they connect SO with ANl (or with AN2) or ANO with SI (or with S2). However, there is a subtle 
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difference: in the case of the branches D1-D4 the couphng of the components ui and U-i did not 
involve any contribution from (or to) component uq, since in this case the last term in the system 
of Eqs. ()2.ip - (|2.2p is inactive. However, in the presence of a nonvanishing component uq, this last 
(four- wave- mixing type) term becomes activated and even though the branch connects, say, SO with 
ANl (as is the case for the branch C2), this waveform introduces a contribution in the component 
U-i, rendering the solution nonzero in all three components. Furthermore, the nature of the profile 
of uq and Ui combined with the functional form of this term (ex UqUi) determine the parity of 
the U-i component [in the case of C2 it should be anti-symmetric as is indeed shown in Fig. [5]. 
Out of this genuinely three-component solution C2, bifurcates an asymmetric variant thereof (again 
through a supercritical pitchfork) destabilizing the solution and the resulting new branch, CI, is 
asymmetric and again involves all three components of the spinor condensate. [It is worthwhile 
to note in passing that as shown in Fig. [5l the solution branch C2 could even be unstable prior 
to the critical point of the symmetry breaking pitchfork bifurcation due to a complex quartet and 
an associated oscillatory instability]. Branches C4 and C3 are the analogs of branches C2 and CI 
but involve predominantly the component U-i (and weakly generate the component Ui). Branches 
C3 and CI eventually merge with AS2 and ASl respectively. Then, branch C6 involves a similar 
coupling which joins S2 and ANO, as is shown in Fig. [T]and Fig. [6] (instead of SO and AN2); out of 
that bifurcates the stable, three-component asymmetric branch C5 also shown in Fig. [6l The analog 
of C6 and C5 involving a coupling between components ui and uq (instead of U-i and uq) is given 
by the branches C8 and C7 which have also been identified in Fig. [T](but for = 0.48). Branches 
C7 and C5 both merge with ASO eventually. Notice that in the rightmost panel of the third row 
of Fig. [H the branch C8 joining SI and ANO and the stable asymmetric branch C7 bifurcating off 
of that again appear to almost coincide (norm- wise). Importantly, it should be noted that none of 
these solutions can be captured in any way by the SMA. 

On the other hand, as is evident in Fig. [U by the comparison of the top right panel of the 
two-mode Galerkin-type approximation with the full results of the original system, the two-mode 
approximation performs very well in providing a semi-analytical prediction (obtained through the 
solution of a few simple nonlinear algebraic equations) of the bifurcation diagram of the full system. 
It is clear that quantitatively the two-mode approximation represents the relevant branches less 
accurately for higher values of the number of atoms N; e.g., the branch CI reaches values of N ^ 3 
as /i ^ 0.52 (at the right end of the Fig. [T|), while from the two- mode approximation it only reaches 
N ^ 2.5. Nevertheless, the method traces in a systematic manner all the branches that can also be 
found in the full numerical results and their corresponding bifurcations. 

From the above, it can be inferred that the full bifurcation diagram is extremely complex in 
the spinor BEC case involving not only "pure" one-component states, but also additional two- 
component and in fact even fully three-component "spinorial" states. The new emergent states 
are always a form of intermediary between original pure states and also possess symmetry break- 
ing bifurcations of their own giving rise to genuinely spinorial, spatially symmetry broken states. 
Despite the fact that many of these states (including the asymmetric ones or even the symmet- 
ric/antisymmetric ones prior to their pitchfork bifurcation points) are genuinely stable dynamical 
states of the F = 1 spinor BEC system, they are not captured in the framework of the SMA. 

Since Fig. [T] was constructed for a particular value of the chemical potential /i-i, it is worth 
commenting on how the relevant features may change upon variation of For this reason. Fig. 
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Figure 5: (Color online) The top left panel shows the maximum real eigenvalue as a function of /io 
for C2 in Fig. [H The other panels show the profiles of wave functions i^i, uq and U-i by solid (blue), 
dashed (green) and dash-dotted (red) lines, respectively, and stability eigenvalues corresponding to 
branches CI and C2 in Fig. [D CI for /io = 0.48 (top right), C2 for /io = 0.453 (bottom left), C2 
for /io = 0.457 (bottom middle) and C2 for /io = 0.47 (bottom right) with /i_i = 0.38. The profiles 
of U-i [dash-dotted (red) lines] are multiplied by 30 to be more visible. 




Figure 6: (Color online) The top left panel shows the maximum real eigenvalue as a function of /io 
for C6 in Fig. [H The other panels show the profiles of wave functions ui, uq^ U-i [by solid (blue), 
dashed (green) and dash-dotted (red) lines, respectively], and stability eigenvalues corresponding 
to branches C5 and C6 in Fig. [JJ C5 for /io = 0.48 (top right), C6 for /io = 0.447 (bottom left), C6 
for /io = 0.45 (bottom middle) and C6 for /io = 0.455 (bottom right) with /i-i = 0.38. The profiles 
of ui [solid (blue) lines] are multiplied by 40 to be more visible. 
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[71 examines both higher and lower values of this chemical potential (both in the full system, as well 
as in the reduced two-mode approximation of the six algebraic equations -left and right panels, 
respectively of the top two rows of the figure-). What is observed is that while a slight increase of 
the parameter does not substantially alter the phenomenology, a slight decrease thereof may result 
in a drastic decrease in the number of branches observed. This can be qualitatively understood as 
follows. When is decreased, initially AS2, later AN2, and finally even S2 will disappear. For 
=0.28 shown in Fig. [7j already branches AS2 and AN2 have disappeared and only S2 persists 
among the pure U-i branches. As a result, all the branches that would "connect" to AN2 are also 
forced to disappear, including, e.g., D2 and Dl, D3, D6 or C4, C3 and C5. Hence, as a rule of 
thumb, decreasing the chemical potentials (in this case of repulsive interactions) generally reduces 
the number of available states that can exist. A guide for potentially identifying such more complex 
combined states is the corresponding presence of "pure" states in the system. Notice that in all 
our considerations herein, we have focused on chemical potentials (i-i < 0.61 such that nonlinear 
modes associated with the second excited state do not become "activated". For higher values of 
the chemical potential, obviously such higher excited states will come into play and the two-mode 
Galerkin-type approximation will no longer be sufficient to address their existence. Nevertheless, in 
considering such higher excited states, we have typically observed them to be unstable and hence 
do not pursue them further herein. 

We should note that we also considered along the same vein the bifurcation diagram for the case of 
the ferromagnetic ^^Rb spinor condensate (recall that in this case 5 = fa/i^s = —4.66 x 10~^ < 0). 
While this difference in the sign of 5 plays a key role in important dynamical phenomena such as the 
modulational instability [33 [151 US] ? nevertheless the structure and main features of the resulting 
state diagram are observed to essentially be the same between the two cases, perhaps due to the 
strong confinement considered here within the realm of the double well potential. 

3.3 Dynamics 

We now provide some representative examples of dynamics of the various branches. We commence 
by considering the branches D2, D4 and D6 respectively in the top, middle and bottom row of Fig. 
[HI We can see that the dynamics and instability of these branches does not strongly couple to the 
0-th component of the spinorial state (and it only does because of the initial small amplitude noise 
seeded in that state, a seeding that takes place in all our simulations in the initial conditions). In 
that sense, the dynamical evolution of these states closely resembles the observations of Ref. [34]. 
In particular, we observe a nearly periodic breaking of the symmetry leading to amplification of 
one of the two components in one of the wells, while typically the other component is amplified in 
the other well [a notable exception to that is the case of D6 where the wavefunctions have the same 
parity, both being anti-symmetric]. 

Roughly similar behavior can be observed in the dynamics of the combined branches involving all 
three spinorial components, such as the ones for C2 and C8 in Fig. [9l as well as the ones for C4 
and C6 in Fig. [TOl In Fig. [9l the solution is predominantly supported in the components ui and 
uq (respectively U-i and uq in Fig. [10]) and has a small amplitude in the remaining component. 
The two predominant components exhibit similar recurrent behavior, whereby at roughly periodic 
intervals the solution becomes asymmetric with stronger support for one component in the one well. 
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Figure 7: (Color online) The norm of the numerically found (left) and the approximate two-mode 
(right) solutions of Eqs. (|2.ip - (|2.2p as a function of /io for M-i = 0.48 (top) and = 0.28 (middle) 
in the case of ^^Na spinor BEC, and = 0.38 in the case of ^^Rb spinor BEC (bottom). The 
notation is the same as in Fig. [H 
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Figure 8: (Color online) Spatio-temporal contour plots of the densities of unstable combined two- 
component solutions in the case of ^^Na spinor BEC. The panels show the simulated evolution of 
wave functions V^i (left), i/jq (middle) and i/j-i (right) in unstable solutions of D2 (top), D4 (middle) 
and D6 (bottom) from Fig. [H respectively. 
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while stronger support for the other component is in the second weh. The weaker component is 
eventuahy more strongly amplified by the instability, but still does not appear to play a critical 
role in affecting the dynamics of the two dominant components. Nevertheless, we clearly see the 
symmetry-breaking manifestation of the relevant instability and the recurrent emergence of the 
ensuing asymmetric waveforms. 



X 01 



0.4 
0.3 
0.2 
0.1 



100 200 300 400 500 
t 



1000 2000 3000 
t 



100 200 300 400 500 
t 



1 


1 


0.1 


lol 






0.08 






0.06 


X ol 






0.04 


1 


|o.02 


lol 



lOl 



0.03 
0.025 
0.02 
0.015 
0.01 

0.005 -10 



X 



0.15 

0.1 

0.05 



1000 2000 3000 
t 



lOl 



-10| 



I 



xlO 
6 



500 
t 



1000 2000 
t 



1000 



X 10 

r 

1.5 
1 

10.5 



3000 



Figure 9: (Color online) Spatio-temporal contour plots of the densities, 



and 



l^-il 



of 



unstable combined three-component solutions in the case of ^^Na spinor BEC. The top and bottom 
panels show the simulated evolution of wave functions ipi (left), tjjo (middle) and tjj-i (right) in 
unstable solutions of C2 (top) and C8 (bottom) from Fig. [H respectively. 
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Figure 10: (Color online) Same as in Fig. [9l The top and bottom panels depict the simulated 
evolution of the densities of components 1 (left), (middle) and —1 (right) in unstable solutions of 
C4 (top) and C6 (bottom) as shown in Fig. [H respectively. 

In all of the above examples, we have seeded the instability by a random (uniformly distributed) 
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perturbation imposed to the original stationary solutions with a small amplitude (10~^), essentially 
emulating the presence of a background of experimental "noise" in a potential experiment. While 
this was not based on a physical model of quantum or thermal noise present in an actual experiment, 
the main thing that matters (at least in as far as the deterministic evolution of the mean field spinor 
models considered herein is concerned) is the projection of the relevant perturbation on the unstable 
eigendirection (i.e., eigenvector) of the linearization around the solution. Our results (for different 
amplitudes, as well as different realizations of the noise) clearly illustrate this fact (see below). 

To illustrate the effect of the noise amplitude, we have repeated our numerical simulations with 
higher noise amplitudes (but with the same spatial distribution of the noisy perturbation); see 
Figs. [11] and [121 We can observe that in such a case the dynamics is principally similar [in all 
cases, the symmetry breaking still occurs, although its evolution in time may differ somewhat for 
different amplitudes cf. FigIT2]. The amplitude of the noise does play a role in the time scale of 
the manifestation of the instability, since the larger the initial noise amplitude, the shorter the time 
interval until it gets amplified (by the instability) to an 0(1) perturbation. 

Finally, to illustrate that the realizations of the noise (for the same noise amplitude) are not sub- 
stantially different, we also did the following numerical experiment. We took ten different (spatial) 
realizations of the randomly distributed perturbation, but all of them with the same amplitude 
and then averaged the result. The relevant findings are presented in Fig. [131 The top panel of 
the figure shows the result of the average of the 10 random realizations of the perturbation (to be 
compared with the corresponding panel of Fig. [9|). The difference between the two panels is shown 
in the bottom of Fig. [131 ^ind is clearly minimal. Hence, for all ten realizations of essentially the 
same noise amplitude, the manifestation of the instability was essentially similar. This validates our 
claim above that the key element in this deterministic dynamical evolution is the projection of the 
perturbation to the dominant unstable eigenmode which dictates the dynamics of the instability 
development and its symmetry breaking manifestations. 

Lastly, although, in the present work, we have used the two mode reduction as merely a tool for 
identifying the static solutions (and their bifurcations) for the full spinor system, let us briefly 
comment in passing about the dynamical properties of the two-mode model. Similarly to what has 
been suggested in [26] , we find (results not shown here) that the dynamical results of the two- mode 
approximation are closer to the ones of the full GP equation for lower nonlinearities (i.e., closer to 
the linear limit) and higher potentials (i.e., a more pronounced double well structure). 

4 Conclusions 

In this work, we considered the statics and dynamics of a F = 1 spinor condensate confined in a 
double well potential. We illustrated that the two-mode Galerkin-type approximation, previously 
developed for one-component [29] and two-component [34J settings can be extended to this 
genuinely three-component setting. The advantage of such a methodology is that it can offer con- 
siderable insight on the full set of stationary states that the system can exhibit, not only as pure 
states involving one-component, or two-component combinations (involving the and ip-i com- 
ponents), but even fully three-component spinorial states. An additional strength of the method is 
that it does not rely on the single mode approximation (SMA) that necessitates the same spatial 
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Figure 11: (Color online) Same as in Fig. [9l The top and bottom panels show the simulated 
evolution of the density of components 1 (left), (middle) and —1 (right) in unstable solutions of 
C2. A random (uniformly distributed) noise of amplitude 10~^ (top) and 10~^ (bottom) is involved 
initially in perturbing the exact unstable solution, compared to the one of amplitude 10~^ in Fig. 
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Figure 12: (Color online) Same as in Fig. [9l The top and bottom panels show the simulated 
evolution of components 1 (left), (middle) and —1 (right) in unstable solutions of C4, with an 
initial random noise of amplitude 10~^ (top) and 10~^ (bottom), compared to the one of 10~^ in 
Fig. M 
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Figure 13: (Color online) Same as in Fig. [H The top panels show the average of 10 simulated 
evolutions of the density of components 1 (left), (middle) and —1 (right) in unstable solutions of 
C2 with different initial noises of amplitude 10~^. The bottom panel is the difference between the 
top panels in this figure and the top panel in Fig. [9l 

profile among all the hyperfine components, but rather it permits to fully explore non-SMA states 
that the system clearly and abundantly possesses. In fact, most of these states are observed to 
be dynamically stable in at least a fraction of parameter space of their existence (i.e., effectively 
for appropriate atom numbers) and hence should be accessible to relevant experiments with spinor 
condensates in double well potentials. We have observed that the two-mode Galerkin-type approx- 
imation is very efficient in unraveling the full bifurcation diagram of the possible states. Finally, 
we have illustrated the dynamics of either two-component or three-component spinorial states in 
direct simulations, observing typically the emergence of the symmetry breaking instability, leading 
to a stronger population of one or the other well, and subsequent recurrence of such asymmetric 
patterns. 

There are many directions that open up with respect to this analytical description of the double 
well system. On the one hand, one can extend this two-mode Galerkin-type approximation to a 
higher- number of mode description (that, for instance, will involve higher excited states) so that 
one can describe the possible situations for higher numbers of atoms/chemical potentials. On the 
other hand, one can extend the present theory to the mean-field description of a F = 2 spinor 
condensate following, e.g., the description of [l^ (see also references therein). Finally, yet another 
interesting direction (even for the F = 1 case) is to extend the present considerations to the 
simplest higher-dimensional case, involving four wells arranged along the nodes of a square and 
attempting a corresponding four-mode reduction for each of the components. In that setting, it 
would be especially interesting to identify multi-pole structures and topological states, such as 
vortices, among others. Such studies are currently in progress and will be reported elsewhere. 
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A Special Case Examples of Solutions 

In this Appendix we consider some special cases, in which the solutions of the algebraic equations of 
our two-mode reduction are analytically tractable (or, in any case, reduce to previously addressed 
problems with a lower number of components). More specifically, there exist one- component states 
which can be found in the case where, e.g., c^^^^^ = 0. In this case, we can readily identify a purely 
symmetric state with 

(ci°))^ = ^^, cr=0, (A.l) 
as well as a purely anti- symmetric state with 



ci°)=0, (4°))^ = ^. (A.2) 



Finally, there is also a mixed or asymmetric state with 



/ (0)n2 _ (MQ - ^a)^b — 3(/iQ - U;b)Tab (0)^2 _ (MQ ~ ^b)^a - 3(/io - C^ajPg^ 



(A.3) 



The mixed branch bifurcates beyond a critical value of /io from the symmetric state if z^^ < 0, 
or from the antisymmetric state if z^s > [29]. It is important to note here that such pure and 
mixed states exist also in components ±1, with the only difference in their definition (except for 
the ^ ±1 in the indices above) that Us ^ + ^a- 

In addition to the above setting, there is another case where these pure and mixed mode states 
appear (now in all three components), namely in the spinor BEC description in the framework of 
SMA [m [371 [38] (see also discussion in [18j). In that case, /io = Mi = M-i = and Uj = Sju{x), 
where u satisfies the single component equation 

jiiu = Cu -\- UsU^ ^ (A-4) 
while the coefficients Sj containing the spin degree of freedom are given by, 

si = cos^(^), 50 = y2cos(^) sin(^), s_i = sin^(^) (A.5) 

si = -^sm{(3), so = cos(^), 5_i = sin(/3), (A.6) 

where is a free parameter. In the context of SMA, which is generally valid if the condensate width 
is far smaller than the spin healing length (so that the terms proportional to Ua do not substantially 
inffuence the dynamics), again the possible solutions reduce to the above pure and mixed modes. 

Finally, yet another special case is the one corresponding to the existence of two-component states. 
In paricular, in the case uq = 0, the three-component system degenerates to the two-component 
setting recently considered in [34], with self-phase modulation proportional to + i^a and cross- 
phase modulation proportional to Us — i^a- Based on the analysis of [34J, two-component states, 
both symmetric and antisymmetric, but also mixed two-component states are expected to exist. 
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